Phase diagram of a Disordered Boson Hubbard Model in Two Dimensions 
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We study the zero-temperature phase transition of a two-dimensional disordered boson Hub- 
bard model. The phase diagram of this model is constructed in terms of the disorder strength and 
the chemical potential. Via quantum Monte Carlo simulations, we find a multicritical line sepa- 
rating the weak-disorder regime, where a random potential is irrelevant, from the strong-disorder 
regime. In the weak-disorder regime, the Mott-insulator-to-superfiuid transition occurs, while, in 
the strong-disorder regime, the Bose-glass-to-superfluid transition occurs. On the multicritical line, 
the insulator-to-superfluid transition has the dynamical critical exponent z = 1.35 ± 0.05 and the 
correlation length critical exponent u — 0.67 ± 0.03, that are different from the values for the tran- 
sitions off the line. We suggest that the proliferation of the particle-hole pairs screens out the weak 
disorder effects. 
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The boson localization due to a random potential 
has continuously attracted significant attention as a 
paradigm of the zero-temperature quantum phase tran- 
sition The superconductor- insulator transition has 
been be ieved to be a manifestation of the boson local- 
ization transition. Two-dimensional realization of this 
transition may be found in disordered thin-film super- 
conductors [0, Josephson-junction arrays 1^, and ''He 
films adsorbed in porous media [Q . An interacting boson 
model, called boson Hubbard model, has been proposed 
1^ to describe the transition between an insulator and a 
superfluid (SF). In the disorder- free case, the insulating 
ground state is a Mott insulator (MI), which has com- 
mensurate boson density and a finite Mott gap which 
suppresses the excitations of freely moving particles or 
holes. As the gap vanishes the system becomes a super- 
fluid. In the presence of disorder, even with vanishing 
energy gap, particles or holes excited can be localized 
by a random potential, resulting in a Bose glass (BG) 
insulator. 

The interplay of interaction and disorder has attracted 
considerable interest in this model. It has been argued 
that in the presence of disorder the transition from the 
insulating to the superfluid state occurs only through the 
BG phase on the assumption that arbitrarily weak dis- 
order is always relevant in two dimensions ||^ . Recently, 
however, quantum Monte Carlo studies of the boson Hub- 
bard model in the Villain representation have shown that 
transitions at or near the tip of the MI lobes, the 
points with the particle-hole symmetry, are most likely 
direct from the MI state to the SF state for weak dis- 
order. Path-integral quantum Monte Carlo simulations 
lH and real-space renormalization calculations of 
the boson Hubbard model support the scenario of the 
direct MI-SF transition at the tips. A simple scaling 
argument combined with renormalization calculation at 
the mean- field level predicts that the direct MI-SF 
is possible around the tip of the MI lobe in high dimen- 



sions {d > 4) while it is possible only at the tip in low 
dimensions (2 < d < 4). Field-theoretical renormaliza- 
tion group studies at the tip, on the other hand, show 
that disorder is always relevant in two dimensions ||l2|,|l3| . 
Strong-expansion studies suggest that the direct transi- 
tion is always unlikely to occur in the presence of dis- 
order I p^ . Similar problems in one dimension have also 
attracted considerable interest recently |l^. Here, it ap- 
pears that a direct MI-SF transition is not supported. 

In this work, we investigate the onset of the superfiu- 
idity in a two-dimensional boson Hubbard model in the 
presence of disorder via quantum Monte Carlo simula- 
tions which employ a (2-|-l)-dimensional classical action. 
We find that a multicritical line exists, separating the 
weak-disorder regime from the strong-disorder regime, 
and on the line, the insulator-to-superfiuid transition is 
associated with novel values of the critical exponents: i.e. 
z = 1.35±0.05 and f = 0.67±0.03, where z and v are the 
dynamical and the correlation length critical exponents, 
respectively. These results are summarized in Figure ^. 
It shows that the direct MI-SF transition survives around 
the tip of the MI lobe for weak disorder. This means that 
the weak disorder is irrelevant for the localization tran- 
sition in two dimensional interacting boson systems that 
have near integer fillings. Strong disorder changes the 
nature of the transition to that of the BG-SF transition. 

We consider a boson Hubbard model on the two- 
dimensional square lattice, given by the Hamiltonian 



(1) 



where bl (bi) denotes the boson creation (destruction) 
operator at site i {rii = bjbi), fii the local chemical po- 
tentials, U the on-site repulsion energy, t the hopping 
strength to the nearest neighbors, and finally the last 
sum is over nearest neighbor pairs. Disorder effect is em- 
bedded in local chemical potential as fii — fj, -\- Vi where 



Vi are random variables independently and uniformly dis- 
tributed in the range [—A, A]. Thus A characterizes the 
strength of disorder. 

To perform Monte Carlo simulations, we follow 
the standard procedure [ p^ to transform the two- 
dimensional Hamiltonian in Eq. (|^) to the (2+1)- 
dimensional classical action 

, V-.7=0 

■^^^ E [J^{x,y,r)-2{p + v,)J^ix,y,T)] , (2) 
{x,y,T) 

where x and y are spatial coordinates and r is imaginary- 
temporal coordinate. Here K ^ ^JijU , analogous to the 
temperature in classical systems, and J{x,y,T) is the 
integer current vectors which measure the fluctuations 
along the corresponding direction of the components. In 
the transformation from Eq. (|l]) to Eq. (^) one assumes 
that superfluidity is destroyed only by phase fluctuations. 

We are interested in the phase diagram of Eq. (||) in the 
space of {K,iJ,,A). The transition described by the clas- 
sical action is studied by Monte Carlo simulations. We 
perform the simulations on the lattice of various sizes, 
denoted hy L x L x L^, where L and Lt are sizes of 
the systems along a spatial and the imaginary-temporal 
axis respectively. The periodic boundary conditions are 
adopted. In order to extract the critical properties of 
the transition, we analyze the data using the finite-size 
scaling theory. An important quantity which indicate the 
onset of the superfluidity is the superfluid stiffness, which 
is measured by the formula 
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where ria; = i'^/L)J2{x,y,T)J''i^^y^'^) the winding 
number along the x direction. Here [...jav denotes the 
average over different realizations of disorder. The finite- 
size scaling behavior of the superfluid stiffness is given 

by [0 



(4) 



where 5 — {K — Kc) is the distance from the critical 
point Kc for each A and /i, p is a scaling function, and 
d is the spatial dimension. Throughout this work, we set 
d — 2. Another useful quantity is the compressibility, 
which directly shows whether the Mott gap vanishes at 
the transition. It is given by the formula 



(5) 



with iIt- — [Ij Lr) y r) ^'^{^1 ''') ^ ^'^'^ assumcd to 
take the scaling form 



with another scaling function k. 



(6) 



In order to investigate the scaling behavior of the su- 
perfluid stiffness and the compressibility, one must spec- 
ify the dynamical critical exponent, z, in advance so as 
to fix the aspect ratio L^/L^ throughout the simulations. 
We have tried various values of z and have chosen the one 
which gives the best scaling behavior satisfying Eq. (jj). 
For lattices with non-integer L^., the superfluid stiffness 
is obtained by a simple interpolation of the two values 
measured in the lattices of nearby integer sizes. 

The effect of disorder at various strength is investigated 
on /i = plane. For strong disorder, say for A > 0.45, 
the onset of superfluidity follows the BG-SF transition 
behavior with z = 2 and v — 0.9, and the compressibility 
is finite at the transition as suggested by the scaling ar- 
gument Jsl and confirmed by subsequent numerical simu- 
lations um. For A < 0.35, however, the transitions take 
the signature of the direct MI-SF transition with z = 1 
and v = 0.67, as reported in the recent simulations 
while the compressibility vanishes at the transition. At 
the intermediate strength of disorder, naive scaling anal- 
ysis using either z = 1 or z = 2 fails. In this work, we 
take a different approach and allow the possibility of in- 
termediate values of z by varying the aspect ratio L-r / 
until the scaling plots collapse onto a single curve. 

A new scaling behavior emerges, at the intermediate 
strength of disorder, which has z = 1.35 ± 0.05 and 
V — 0.67 ± 0.03. This scaling behavior is shown in Fig- 
ure H We observe the best scaling when z = 1.35. In ad- 
dition, with this value of z, KcH obtained from the scal- 
ing behavior of the superfluid stiffness and the compress- 
ibility data curves are consistent with each other. The 
compressibility vanishes at the transition as expected for 
z < d (see Figure ^(b)). Kisker and Rieger reported 
Q z « 1.4 for the intermediate strength of disorder, 
A = 0.4. However, they interpreted the result as sim- 
ply measuring an effective exponent. As we can see in 
Figure ^ the robust finite-size scaling behavior for vari- 
ous sizes strongly suggests that the transition truly has 
new values of the critical exponents. In order to rule 
out the possibility that we are disguised by finite size ef- 
fects, we have performed simulations adopting a different 
aspect ratio and find the same critical exponents. The 
value of the correlation length critical exponent happens 
to be the same as the pure transition at the tip. 

This behavior persists even off the tip (/i ^ 0). The off- 
tip transition of the pure case, usually called the generic 
transition, is of a mean-field type, so that z = 2 and 
1/ = 1/2 for the MI-SF transition. If disorder is rele- 
vant, we expect in general that the correlation length 
critical exponent have a different value |18|. Thus, we 
can numerically identify the BG-SF transition from the 
MI-SF transition by measuring ly as well as by observ- 
ing the compressibility. Previously Park et al. have re- 
ported that even in the presence of weak disorder, the 
generic transition survives near the tip. As A increases, 
the generic MI-SF transition ends at the critical value, 
Ac(^), above which the BG-SF transition occurs. 

The critical strength of disorder, Ac(/x), therefore, de- 
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fines a multicritical line, separating the critical surface 
KciUjA) into the two regions: the strong-disorder re- 
gion and the weak-disorder one. Figure ^ is the measured 
Ac(/i). Numerically we find good scaling behaviors us- 
ing the above exponents for a range of parameters (see 
the inset of Figure ||). The error bars in Figure || show 
such ranges. We believe that, in larger systems, the finite 
range should shrink to a point. Otherwise, the existence 
of an incompressible localized phase other than the BG 
phase is required, which is unreasonable. 

The existence of the multicriticality strongly suggests 
that in the presence of disorder two different universal- 
ity classes exist, in one of which disorder is irrelevant. 
The irrelevance of weak disorder may be understood from 
the abundance of the particle-hole pairs excited. At the 
tips, the proliferation of particle-hole pairs near the tran- 
sition modifies the particle propagation to yield z = 1 
(otherwise z = 2) in the pure case. These particle-hole 
pairs may screen the disorder effects. In order to directly 
check the possibility that the proliferation of particle- 
hole pairs may screen out the weak disorder effects, we 
have repeated Monte Carlo simulations of the classical 
model given by Eq. (^ under the condition that the 
background particle-hole fluctuations are completely sup- 
pressed. For simulations of the system governed by the 
action of Eq. (|^), one should generate all possible current 
configurations {J(x,y,r)} |16|. The suppression of the 
particle-hole fluctuations is realized by disallowing any 
Monte Carlo updating which generates local particle-hole 
current loops. In the spatial planes, however, all possi- 
ble current configurations are allowed as usual. Thus the 
allowed current configurations, if they are not restricted 
in a spatial plane, percolate in the r-direction (the imag- 
inary temporal direction). 

Figure ^ shows the scaling behavior of the super- 
fluid stiffness when the particle-hole fluctuations are sup- 
pressed at the tip {fj, = 0) for relatively weak disorder 
(A = 0.3). The scaling shows that the critical exponents 
are z = 2 and ly = 0.9, indicating the BG-SF transition 
occurs. Note that, with the particle-hole fluctuations, the 
transition has the behavior of the MI-SF transition [^|J^] 
at this strength of disorder. This simulation result, there- 
fore, are consistent with the argument that the particle- 
hole pair fluctuations make the weak disorder irrelevant 
to the transition. 

Now one may question how the particle-hole fiuctu- 
ations screen the disorder. In the mean-field picture, 
ignoring the particle- hole fiuctuations, particles (holes) 
are localized by arbitrarily weak disorder due to the co- 
herent impurity back-scattering. However, the prolifera- 
tion of the particle-hole pairs could hinder the coherent 
back-scattering. A particle propagating could be cap- 
tured and destroyed by holes, and a new particle starts 
to propagate. A phase shift which accompanies such 
particle-exchange process may destroy the coherent back- 
scattering of the particle. As a result, this effect could 
screen out a weak random potential, making it irrelevant. 



Slightly off the tip, there are still abundant particle-hole 
pairs fluctuating. Thus this screening effect will survive 
even off the tip. 

In summary, via quantum Monte Carlo simulations of 
the phase-only model of a disordered boson Hubbard 
model, we have found that there is a multicritical line 
which divides the universality of this model into the 
weak-disorder regime and the strong-disorder regime. In 
the weak-disorder regime, disorder is irrelevant for the 
insulator-to-superfiuid transition. On the multicritical 
line, the dynamical critical exponent z ~ 1.35 ± 0.05 and 
the correlation exponent v — 0.67 ± 0.03 are obtained. 
Our result strongly suggests the direct MI-SF transition 
survives for weak disorder. We argue that the particle- 
hole pairs fluctuating near the tip will screen out disorder 
to make a weak random potential irrelevant. 
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FIG. 1. Schematic phase diagram of a disordered boson 
Hubbard model. Here jj, is the chemical potential, A is the 
strength of disorder, and K is the parameter characterizing 
the hopping of bosons. The particle-hole symmetry is defined 
by the condition ^ — 0. The multicritical line separates the 
critical surface into the direct MI-SF transition region and 
the BG-SF transition region. 
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FIG. 2. The finite size scaling behavior of (a) the super- 
fluid stiffness and (b) the compressibility at the tip of the MI 
lobe {fi = 0) with A = 0.40. We set z = 1.35. The curves 
cross at Kc ~ 0.292. The inset shows the scaled data along 
A'-direction with u — 0.67. 




4 



FIG. 3. The critical strength of disorder, Ae(/Lt), as a func- 
tion of the chemical potential /i. At the tip {n = 0), we have 
tuned the value of A to find the critical one. The error bar 
indicates the range in which the scaling plots show the be- 
havior expected on the multicritical point. Ofi' the tip fi is 
tuned instead, while A is fixed. Inset: The superfluid stiff- 
ness scaling curves away from the tip (yu = 0.175) at A = 0.2 
with z = 1.35 and v = 0.67. 
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FIG. 4. The finite size scaling behavior of the superfluid 
stiffness at A = 0.3 {jj, = 0) for the systems in which the 
particle-hole pair fluctuations are completely suppressed. We 
set z = 2. The inset shows the scaled data along /C-direction 
with v = 0.9. These exponents indicate that the nature of 
the transition is of the BG-SF one. 
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